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ABSTRACT 

We have employed a 3-d energy-conserving radiation MHD code to simulate 
the vertical structure and thermodynamics of a shearing box whose parameters 
were chosen so that the radiation and gas pressures would be comparable. The 
upper layers of this disk segment are magnetically-dominated, creating condi- 
tions appropriate for both photon bubble and Parker instabilities. We find little 
evidence for photon bubbles, even though the simulation has enough spatial reso- 
lution to see them and their predicted growth rates are high. On the other hand, 
there is strong evidence for Parker instabilities, and they appear to dominate the 
evolution of the magnetically supported surface layers. The disk photosphere is 
complex, with large density inhomogeneities at both the scattering and effective 
(thermalization) photospheres of the evolving horizontally-averaged structure. 
Both the dominant magnetic support and the inhomogeneities are likely to have 
strong effects on the spectrum and polarization of thermal photons emerging 
from the disk atmosphere. The inhomogeneities are also large enough to affect 
models of reflection spectra from the atmospheres of accretion disks. 

Subject headings: accretion, accretion disks — instabilities — MHD — radiative 
transfer 
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Introduction 



The first attempts to understand accretion disks concentrated on deriving radial pro- 
files of a zimuthally-averaged a nd ve rtical ly-integrated quanti t ies in time-steady disks. For 
example, IShakura Sunyaevl (119731 ) and iNovikov fc Thornd (119731 ) worked out the theory 
for such disks around blac k holes by applying the con straints of energy and angular momen- 
tum conservation. When IShakura fc Sunyaevl (119731 ) introduced the notion that the stress 
responsible for moving matter inward through a disk should be proportional to the pressure, 
it was with a view toward tying the stress to local conditions in an azimuthally-averaged and 
vertically-integrated sense. Efforts were made within that framework to compute the actual 
vertical thickness as a function of radius, but it was recognized that this computation was 
subject to considerable uncertainty due to lack of knowledge of internal disk dynamics. 

Particularly with the recognition that disk stresses can largely be explained by correlated 
M HD turbulence stirred b y the magneto-rotational instability (as reviewed, for example, 



in 



Balbus fc Hawleyl Il998l ). recently attention has been given to what may happen inside 



disks, even when their vertical scale-heights are quite small compared to the distance to 
the central object. This attenti on has included both analytic work on wave modes and 
instabilities that may exist (e.g. iNowak fc Wagoner! Il993l ; iGammid 1 19981 ; iKim fc Ostriker 
2 OOOl ; iPessah fc Psaltia I2005T) and numerical simulations of vertically stratified shearing-box 
segments of disks ([Brandenburg et al.l Il995l ; IStone et al.l Il996l ; iMiller fc Stond |2000| ; iTurner 
20041 : iHirose. Krolik. fc Stone! l2006h . 



Special problems arise in understanding the vertical structure of disks when radiation 
pressu re is comparable to or greater than gas pressure. As first pointed out by lShakura fc Sunyaev 
( 119731 ). this condition is likely to occur in the inner portions of disks accreting at more than 
a small fraction of Eddington, particularly when the central mass is 3> 1M . Perhaps the 
biggest potential problem is that thermal and viscous instabilities in th e vertically-integrated 
quan t ities might occur as a re sult of radiation pressure dominance (ILightman fc Eardley 
1974 ; IShakura fc Sunyaevl Il976l ). 



The large radi ation pressure regime may also be vulnerable to smaller scale photon 
bubble instabilities (lGammidll998l ). The most rapidly-growing form of this instability occurs 
at short wavelengths, where photons diffuse rapidly, resulting in radiative amplification of 
magnetosonic waves (IBlaes fc Socrates! 120031 ). At least for strong magnetic fields, the initial 
nonlinear evolut i on of the instability in this shor t wavelength limit produces trains of shocks 
(lBegelmanll200ll ; ITurner et al.ll2005l ). iBegelmanl (120021 . 120061 ) has suggested that these shock 
trains would permit accreting black holes to radiate at luminosities 1-2 orders of magnitude 
larger than Eddington. 
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We are embarked on a program to explore the nature of the internal dynamics in disks at 
varying levels of relative importance of radiation pressure support. We do this by conducting 
detailed numerical simulations that both conserve energy to high accuracy and include radia- 
tio n transport in a vertically stratified disk annulus. The first such simulation was performed 



by iTurnerl (120041 ) with a high ratio of radiation to gas pressure at the midplane. While very 
interesting results were obtained, the simulation achieved only partial energy conservation 
and was not able to reach low enough densities to incorporate the photosphere within the 
simulation domain. Solutions to these problems have since been found and a simulation 
of a disk an nulus in which the radiation p ressure was about 20% of the gas pressure was 



prese nted by iHirose. Krolik. fc Stone! (120061 ) . In a companion paper to this one (IKrolik et al. 



20071 ). we report on the thermodynamic properties of a simulated disk annulus in which the 



time-averaged radiation and gas pressures were roughly equal. 

In this paper, we investigate internal structures near the surface of the same disk annulus 
whose thermodynamics were studied in the companion paper. We concentrate on these upper 
regions because the shorter length scale instabilities predicted by linear theory tend to have 
their fastest growth rates there. As we shall show, two are of special interest: the photon 
bubble instability and the Parker instability. In addition, it is the structure of these surface 
layers that is most important in determining the degree to which the photon spectrum 
emerging from the disk deviates from a black body at the local effective temperature. 

This paper is organized as follows. In section 2 we briefly summarize the methods 
and parameters of the simulation we are investigating. Then in section 3 we describe the 
properties of photon bubble and Parker instabilities that we expect based on applying linear 
instability theory to the horizontally averaged structure in the simulation. In section 4 we 
compare these expectations to what we actually observe in the three-dimensional structure of 
the simulation. Strong inhomogeneities are observed near the photosphere in the simulation, 
both in density and, to a lesser extent, temperature, and we describe these in section 5. In 
section 6, we discuss the possible implications of our results for models of the emergent and 
reflected photon spectra, and for the emergent photon polarization. Finally, we summarize 
our results in section 7. 



Calculational Method 



Our results are based on a 3-d radiation MHD simulation of a vertically stratified section 
of an accretion disk, computed under the shea r ing-b ox approximation. The code is described 
in greater detail in IHirose. Krolik. fc Stond (120061 ). although the version we used differs 
slightly from that one by certain technical improvements. In its essential elements, this code 
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computes the equations of 3-d ideal MHD, capturing all grid-scale numerical dissipation as 
heat. Radiation transport and radiation forces are included, but with thermally- averaged 
opacity and emissivity and in the flux-limited diffusion approximation. Further details about 



initia l conditions and boundary conditions are given in the companion paper (IKrolik et al. 



20071 ). which concentrates on describing the thermodynamics and larger-scale properties of 
this simulation. 

In order to arrive at a disk annulus in which the radiation and gas pressures are compa- 
rable, we chose the following parameters for this shearing-box segment: The central mass is 
6.62M , the radius is r = 150r s = 1.5 x 10 8 cm, and the surface density is 4.7 x 10 4 g cm -2 . 
If it were described by a Shakura-Sunyaev stress ratio a = 0.02, the mass accretion rate 
would produce a luminosity O.lLg if the efficiency were 0.1. The local effective tempera- 
ture would then be 9.0 x 10 5 K. Our unit of l ength H = 3 .1 x 10 6 cm is the (half)-disk 



thickness in the gas-dominated limit, as given in iKroliki (119991 ). The midplane optical depth 
is (1.58 ± 0.01) x 10 4 , varying slightly due to the temperature-dependence of the free-free 
contribution. At all times, the opacity is dominated by Thomson scattering. The initial 
magnetic field is a twisted azimuthal flux tube with no net poloidal flux and energy density 
1/25 the midplane gas plus radiation pressure. 

Our computational box had a vertical thickness of 12H, stretching symmetrically above 
and below the midplane. Its radial thickness was (3/4)H and its azimuthal thickness was 
3H. These dimensions were described by 512 x 32 x 64 cells, respectively. The simulation 
ran for 318 orbits, equiva lent to more than 40 thermal times. As discussed in the companion 
paper (IKrolik et al.ll2007l ). the epoch with highest total energy content was at 90 orbits, and 
the epoch with lowest total energy content was at 150 orbits. We focus our analysis in this 
paper on just these two epochs. 



3. Expectations Based on Linear Instability Theory 

3.1. Photon bubbles 

To check for the presence of photon bubble instabilities in the present simulation, we 
must first delineate the regions where photons diffuse rapid ly relative to the tim e scales 
relevant for acoustic waves in a gas and radiation mixture (jAgol fc Krolikl Il998l ). When 
gas and radiation exchange heat rapidly, an excellent approximation for most, but not all, 
regimes of plasma acoustic behavior in the case under consideration here, it is convenient to 
write the dispersion relation for compressive waves in terms of an effective sound speed (e.g. 
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Appendix B of iBlaes fc Socrateal2003l ) 



u [e + AE] cl + ^-AEc 



» t e + 4 ^ + iS~A E 



3kfp 



Here k = 2ir/\ is the wavenumber of the disturbance, uj = kC s is the corresponding angular 
frequency, e is the internal energy density of the gas, E is the radiation energy density, p is 
the density, and Kp is the flux mean opacity. The quantity c x = (pg/p) 1 ^ 2 is the isothermal 
sound speed in the gas, where p g is the gas pressure. The total adiabatic sound speed c t of 
the fluid is defined by 



n2 _ 1QE 2 + 60( 7 ~ l)#e + 97(7 ~ l)e 2 _ Ti{p s + E/3) 
Ct ~ 9(e + AE)p 



P 



(2) 



where 7 = 5/3 is the adiabatic index of the gas and Ti is the first generalized adiabatic expo- 
nent rela ting pressure and volu me changes when radiation and matter are tightly thermally 
coupled (jChandrasekharl 119671 ) . Equation (CEJ) shows that the wave speed ranges from q to 
c t ; the frequency for a given wavelength may therefore be easily bounded. 

It is then straightforward to show that acoustic disturbances with wavelengths satisfying 

2 
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= Ar 



(3) 



are in a regime where photons diffuse rapidly compared to the wave period. As a result, the 
radiation pressure response is lost, and the sound speed approaches the isothermal sound 
speed in the gas c Y . On the other hand, acoustic disturbances with longer wavelengths 
satisfying 



A > A s 



R 



r x 1 + 



E 

3p g 



1/2 



A 



R 



(4) 



trap photons, and therefore propagate at c t , the total sound speed in the gas and radiation. 
For intermediate wavelengths with Ar < A < As, the damping of the disturbances is maximal. 
However, when the radiation energy density is not much larger than the gas pressure, as in 
the present simulation, the maximal damping rate is less than the mode period, and c t ~ q. 

Fig. [T] shows the critical wavelengths Ar and As as a function of height at the highest (90 
orbits) and lowest (150 orbits) total energy epochs in the simulation, based on horizontally 
averaged data. Also indicated is the range of wavelengths in the azimuthal direction that 
can be resolved by the code, from eight azimuthal cell sizes to the box size in the azimuthal 
direction. Near the midplane, photons are trapped at all epochs for acoustic disturbances 
at all wavelengths resolved by the code. Outside ±2.5H, compressive disturbances are in 
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the rapid diffusion limit for all accessible wavelengths at t — 150 orbits, when the disk is 
relatively cool. Even at t = 90 orbits, when the disk is relatively hot, there is a factor of 
two range of resolved wavelengths that are in the rapid diffusion limit in the subphotosphere 
surface layers. 

In addition to the rapidity of photon diffusion, another issue of rel evance to photon 



bubbl es is how rapidly gas and radiation exchange heat. As discussed in iBlaes fc Socrates 



(120031 ) . the gas and radiation temperatures are effectively locked together in a rapidly grow- 
ing photon bubble disturbance provided uj^h = AEKpc/e is greater than the gas acoustic 
frequency g/ci associated with a gas pressure scale height. (Here k is the Planck mean opac- 
ity associated with the free- free absorption used in the code.) The ratio of these two rates 
is shown in Fig. [2J and except for the uppermost layers at 150 orbits, rapid heat exchange is 
a good approximation for photon bubbles. 

In this thermally locked, rapid diffusion regime, photon bubbles reach their maximum 
growth rates 7^ over a finite range of wavelengths. The approximate lower bound is a short 
thermal "cutoff wavelength" A cuto fr = 2-kc-J (cothld) 1 ^ 2 below which rapid heat exchange fails 
to be a good approximation in the amplifying acoustic wave, and the wave stabilizes due 
to the damping associated with the slow heat exchange. The approximate upper bound 
is a long "turnover wavelength" A tU rnover = 7d/^ P h above which the wave is no longer a 
simple magnetosonic wave propagating with phase speed t> p h, and the growth rate declines. 
These wavelength bounds are shown in Fig. [3] for the horizontally-averaged conditions in the 
simulation. 

Fig. H] shows the expected m aximum growth rate 7^ of the photon bubble instability 



(Eq. 93 of IBlaes fc Socrates! 120031 ) as a function of height in the simulation box at 90 and 
150 orbits. We assumed that the relevant wavelengths are in the rapid diffusion limit, which 
Figs. [T] and [3] show to be valid at least in the subphotospheric surface layers. We also 
assumed a purely azimuthal magnetic field and a wave vector orientation inclined upward at 
45 degrees to the horizontal in the vertical/ azimuthal plane. Much of the field is in fact in the 
azimuthal direction in the upper layers and 45 degrees is close to the inclination that gives 
maximal growth. Apart from particular cases where the growth rate strictly vanishes, other 
field orientations and wave vector directions alter the growth rates by geometric factors of 
order unity. It is important to note that we calculate growth rates only inside the Rosseland 
mean photosphere of the horizontally averaged structure, as the photon bubble instability 
does not exist when the medium is optically thin. 

Figs. 3 and 4 show that at t = 90 orbits, modes with wavelengths shorter than about a 
scale height should grow by many e-foldings per orbital period in the surface layers beneath 
the Rosseland mean photosphere. Our simulation has plenty of resolution to see these rapid 
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growth rates: a scale height projected through 45° corre sponds to 15 zones in the ^-direction 
and 30 zones in the ^-direction. Resolution studies by iTurner et al.l (120051 ) found that the 
growth rates predicted by linear theory should be recovered in a numerical simulation to 
within ten percent provided there are more than 15 zones per wavelength. On the other 
hand, Fig. 3 also shows that wavelengths more t han a factor of two sm aller than the 
tu rnover wavelength a re not well-resolved. Figs. 6 of iBlaes fc Socrates! (120031 ) and 12 and 13 
of ITurner et al.l (120051 ) show that the growth rate at the turnover wavelength can be smaller 
by as much as a factor of two compared to the maximum growth rate 7^, but that still gives 
4 — 5 e-foldings per orbital period. The maximum growth rates are considerably smaller 
at t — 150 orbits, and will be reduced further above z = 2.5H in the upper half of the 
simulation due to the breakdown in good thermal coupling between the gas and radiation. 
Photon bubbles are therefore most likely to be present at t = 90 orbits. 



3.2. The Parker instability 



All analyses of the magnetoson ic photon bubb l e instability so fa r have assumed a uniform 
magnetic field in the equilibrium (lGammielll998l ; iBegelmanl l200ll ; IBlaes &: Socrates! 120031 ). 
Even if that field is strong, it does not support the ba ckground against g ravity. In contrast, 
we have shown in section §3.2 of the companion paper (jKrolik et al. 2007 ) that the magnetic 
field in our simulation provides the dominant support against gravity at h igh altitudes. This 
fact immediately suggests that the Parker instability (|Parkerlll966l . 119671 ) might play a role 
instead of or in addition to the photon bubble instability. 

We can get a rough idea of the relevance of the Parker instability to this simulation by 
examining its dispersion relation in an isothermal atmosphere. We do this in the Appendix. 
Figs. [5] and [6] show the maximum growth rate [equation (IA16I) ] and characteristic wavelengths 
[equations (IA14l) - (IA15j) ] of the Parker instability as a function of height at t = 90 and 150 
orbits. These figures are based on horizontally averaged data at each epoch, and on the 
assumption that the magnetic field is azimuthal. We also assumed that radiative diffusion is 
rapid, so that perturbations in the gas are isothermal, though we believe that the results are 
approximately valid even under conditions where the radiation is trapped. This is because 
the Parker instability arises from a competition between magnetic pressure and thermal 
pressure, and the gas and radiation pressures are comparable here. Fig. E] leads us to expect 
that the regions outside approximately ±2H will be Parker unstable during the high total 
energy epoch at t = 90 orbits. In the low total energy epoch at t — 150 orbits, we expect 
the regions z < —2H and H < z < 3H to be Parker unstable. 



The Parker instability growth rates exceed the orbital frequency by roughly a factor of 
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two, and are larger than the maximal photon bubble growth rates except perhaps very near 
the Rosseland mean photosphere at t = 90 orbits. Interestingly, the wavelength for maximal 
growth of the Parker instability roughly corresponds to the azimuthal size of the box, and 
our analysis in the Appendix indicates that the corresponding wave vector should be in the 
horizontal direction. 



4. To what extent are instabilities present in the simulation? 

4.1. Photon bubbles 

The simulation results are somewhat ambiguous as to the presence of the photon bubble 
instability. Fig. [7] shows the density, magnetic field, and velocity field in the vicinity of the 
upper photosphere at 90 orbits, the time of highest total energy in the simulation. Two 
photospheres are actually indicated in this figure, both calculated from the horizontally 
averaged structure at this epoch: the Rosseland mean photosphere and the effective, or 
thermalization photosphere. The latter is the outermost surface where photons exchange 
energy with matter. It is therefore the surface inside of which photons can come into thermal 
equilibrium with the matter. Because Thomson scattering dominates the free-free absorption 
opacity throughout the domain of the simulation, the effective photosphere is significantly 
deeper than the Rosseland mean photosphere. 

In the region 3H < z < AH there are 2-3 relatively high density regions which are 
somewhat reminiscent of the s hock trains that form in the initial nonlinear development of 



the photon bubble instability (iBegelmanl 1200 ll ; Turner et al.ll2005l ). The high density streak 



stretching diagonally from the upper left (— 0.5H,3.8H) to the lower right (0.7H,3.2H) in 
the center of the figure is almost certainly a shock, as the velocity field has significant 
convergence in this region. Two other shock-like features in the lower left of the figure 
[ranging from (-0.9H ,3.2H) to (~0.5H,3.0H) and from (-1.1H,3.1H) to (-0.7F,2.8#)] 
run parallel to the first shock, although whether these are true shocks is not completely clear 
as the velocity field does not show as much convergence as at the first shock. The magnetic 
field lines show a fairly discontinuous change in orientation across the first shock, which 
might be indicative of some buckling under the weight of the dense material. There is also 
a rapid change in orientation of the magnetic field lines near the other two putative shocks. 
Examination of nearby radial slices at the same time shows that the first shock extends 
over only ~ 0.2H in radius. Our simulation data dumps lack the time resolution to explore 
its evolution: it is absent and the overall structure is quite different one orbit earlier and 
one orbit later. This is not surprising, however, as the background differential rotation will 
shear out a region of radial thickness x = 0.2H to extend over 3ixx ~ 2H in the azimuthal 
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direction after one orbit. 

It is unclear whether these shocks are the result of photon bubble instability. To try 
to resolve this ambiguity, we have investigated how the radiation flux and fluid velocity are 
correlated in the layers where the photon bubble instability is expected. At least in the rapid 
diffusion limit of interest here, photon bubbles are fundamentally a radiative amplification 
of a magnetosonic wave. The density compressions and rarefactions in the wave alter the 
diffusion paths of photons in such a way as to produce radiation flu x perturbations tha t 



everywhere make acute angles with the fluid velocities in the wave ([Turner et al.l 120051 ) . 
Oscillating radiation pressure forces in phase with the oscillating fluid velocities are the 
result, which drive the wave to larger amplitudes. To check for this distinctive directional 
relationship between the perturbed forces and velocities, we show in Fig. [S] the distribution 
of cosines of the angle 9 between the fluctuating radiation flux and fluid velocity vectors at 
all grid points in a set of horizontal slices through the simulation at 90 orbits. To be precise, 

_ (F- < F z > z) ■ (v + f fisy) 
C ° SU ~ |F-<F 2 >z||v + fn*y| ' (5j 

where the angle brackets refer to a horizontal average at the particular height in question. 
We subtract off the horizontal average of the vertical flux component because a radiation 
driven outflow would also give rise to a correlation between flux and velocity without having 
anything to do with photon bubbles. 

As can be seen in Fig. [HI at three of the four sample heights, there is essentially no 
directional correlation between the fluctuating part of the radiation flux and the fluid velocity. 
Only at z = +2AH is there a significant positive correlation; in fact, there is a weak anti- 
correlation at z = +3.3H and 3.8H. Ironically, our application of the linear theory of photon 
bubbles to these data suggests that they grow only for z > 2.5H. We conclude, therefore, 
that there is little evidence for photon bubble-driven dynamics here. 



4.2. Parker instability 

Although evidence for photon bubbles is scant, evidence for Parker instability is strong. 
This mode involves an interplay between support against gravity provided by magnetic 
pressure forces and restraint against magnetic buoyancy provided by magnetic tension forces. 
We therefore begin this part of our analysis by mapping out the strengths of these forces. 

The vertical accelerations produced by magnetic tension and magnetic pressure are 

1 / dB z dB z \ 
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and 

*P = ^(Bl+B 2 y ), (T) 
8npoz y 

respectively. These accelerations are plotted in Fig. [9j In a volume-averaged sense, the 
magnetic pressure force almost always pushes vertically outward; by contrast, in the same 
sense of averaging, the magnetic tension force generally holds matter in. In this volume- 
averaged sense, the two magnetic forces appear almost to balance each other. As we shall 
see shortly, in Parker modes there is a strong correlation between the sense of curvature of 
the magnetic field (which controls the direction of the tension force) and the gas density: 
relatively low density and inward tension force are often found together. For fixed total mass, 
low density regions must always occupy more total volume than those of high density, so 
inward magnetic tension regions dominate in a volume average. On the other hand, support 
against gravity is best thought of in terms of mass-weighted quantities (right-hand panel of 
Fig. [9]). Averaged in this fashion, magnetic tension becomes much less important, and we see 
that the upper layers of the disk are, as previously noted, primarily supported by magnetic 
pressure gradients. 

The fact that large inward magnetic tension forces begin outside z = ±3H is an indica- 
tion that the Parker instability is in operation at those altitudes. This height is not far from 
the z = ±2H boundaries where we expected the Parker instability to develop on the basis 
of our simplified linear analysis. 

A further diagnostic of Parker instability can be seen in the density and magnetic field 
structure shown in Fig. [7] at t = 90 orbits. The densest regions of gas lie in troughs of 
magnetic field lines, whereas the lighter regions are located near magnetic field line crests. 
To quantitatively check this visual impression, the left hand panel of Fig. [TU] shows the 
magnetic tension force per unit volume versus the density at each grid zone in a z = 2.75H 
horizontal slice through the box at this epoch. There is a clear correlation between these 
two quantities: gas heavier (lighter) than average is predominantly associated with upward 
(downward) tension forces. 

Fig. [11] shows the upper surface layers at t — 150 orbits. The structure above the 
photosphere is very complex with large flow velocities, and the 4> — z projection of the field 
shows circulation, including regions of nearly vertical field lines. This complex structure is 
no doubt related to the l arge deviations fro m hydrostatic balance (see right panel of Fig. 6 



in the companion paper, iKrolik et al.l 120071 ). Nevertheless, a long wavelength Parker mode 
appears to be present in the region 2H < z < 3H, which is again where it was expected 
to occur based on the minimum unstable wavelength shown in Fig. [6] Its wavelength and 
horizontal orientation also conform with the linear analysis expectations. The right hand 
panel of Fig. [TUJ shows that tension is again highly correlated with density in the z = 2H 
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horizontal slice. 

This correlation between density and tension forces is perhaps the strongest evidence 
that Parker modes play a significant role in the outer layers of the simulation. A volumetric 
rendering of the correlation is shown in Fig. [12] Regions that are denser than average at 
a given height are mostly associated with field line valleys, whereas low density regions are 
mostly associated with crests in the field lines. 

Fig. [T3] shows the fractional horizontal surface area at each height for which the corre- 
lation between density and tension forces holds. Within approximately two scale heights on 
either side of the midplane, there is no correlation. Outward and inward tension forces in 
this region are equally likely to be associated with higher and lower than average densities, so 
the fractional surface area exhibiting the expected behavior of the Parker instability is fifty 
percent. Further away from the midplane at t = 90 orbits, however, the vertically outward 
(inward) tension forces are more likely to be associated with lower (higher) than average 
density, as expected from the Parker instability. This vertical dependence of the magnetic 
tension-density correlation is completely consistent with our linear instability expectations 
from section 3.2 above: unstable growth at \z\ ^ 2H and no instability closer to the mid- 
plane. The situation is approximately the same at t = 150 orbits, except that there is no 
evidence of the expected Parker correlation below z ~ —3H, or near z = 3H, AH, and above 
5H. Our linear analysis of section 3.2 predicted that at this epoch the regions z < —2H 
and H < z < 3H would be Parker unstable. Again, we suspect that the discrepancy arises 
partly from the breakdown of hydrostatic equilibrium in the outer layers at this epoch, as 
shown in the right panel of Fig. 6 of the companion paper. 

It is fortunate that our simulation domain just happened to be the right size for the 
maximal Parker instability wavelengths to fit inside the box. Whether these modes can 
fit in the computational domain is clearly something to be borne in mind when choosing 
parameters for future stratified shearing box simulations. It is unclear whether the finite 
azimuthal extent of the box might be affecting the nonlinear saturated state of the instability. 
One could imagine, for example, that the tension forces could be reduced if the nonlinear 
development of the instability cascaded toward longer wavelengths, but this would not be 
possible in our finite shearing box. Clearly further work with larger simulation domains is 
warranted in the future. 

For comparison, the upper su rface layers in the middle of the box near the end of the gas 



pressure dominated simulation of iHirose. Krolik. &: Stone! (120061 ) are shown in Fig. O The 



structure in the lower regions of this figure again appears to be consistent with the Parker 
instability. Fig. [15] shows that the expected correlation between density and tension forces 
is again widespread outside approximately two scale heights away from the midplane. Not 
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surprisingly, the Parker instability appears to be generically important in the magnetically 
supported upper layers of accretion disks even when there is little radiation pressure. 

5. Photospheric Irregularity 

Very large, irregular density inhomogeneities are evident in the surface layers and extend 
as deep as the effective photosphere of the horizontally-averaged structure at all epochs. The 
distribution of densities at the upper effective and Rosseland mean photospheres at 90 and 
150 orbits is shown in Fig. [16] The density fluctuations range over factors of approximately 
ten above and below the mean, with small isolated rarefied regions that are even less dense 
still. The vertical distribution of horizontally-averaged root mean square fractional fluctua- 
tion of density is shown in Fig. [17] for the same epochs. The large density inhomogeneities 
clearly extend over a broad range of heights outside ±2H. 

One possible explanation for the presence of these density inhomogeneities lies in the 
fact that, in isolated parts of the most rarefied regions outside the effective photosphere, the 
gas temperature can greatly exceed (by greater than a factor 100) the effective temperature 
of the radiation, as shown in Figs. [18] and [19] The code assumes that gas and radiation 
exchange heat purely through free-free absorption and emission, and the density is so low 
in these regions that the gas and radiation are unable to equilibrate through this channel. 
Including Compton scattering in the gas and radiation energy equations would enhance the 
thermal coupling. The ratio of time scales of Compton to free-free cooling of the gas is 
approximately 

^Comp ^ K (^gas — ^~rad) m eC 
tff 4fc B (^gas — T za d)KTT zgd ' 

where k is the Planck mean free-free opacity, kt is the Thomson opacity, T gas is the gas 
temperature, and T ra d is the effective temperature of the radiation. In the conditions of 
the current simulation, free-free is faster than Compton in exchanging heat between the gas 
and radiation throughout most of the volume. However, in the isolated regions where the 
gas temperature significantly exceeds the radiation temperature, the Compton cooling time 
can be shorter than the free-free cooling time by factors of 5 to 10. The gas may therefore 
be artificially too hot in these regions, because we did not include Compton cooling in the 
simulation. Artificially elevated temperature might enhance the degree of inhomogeneity 
by artificially raising the gas pressure. Although this effect may explain the regions we see 
where the gas and radiation temperatures are far apart, it cannot explain all regions of large 
density inhomogeneity: Figs. [TS] and [TH] show that the radiation and gas temperatures are 
equal at the effective photosphere (as they must be), and yet there can be large density 
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contrasts there. 



Photon bubbles would produce inhomogeneities, but we have failed to find any other 
evidence for these instabilities in the simulation. Moreover, substantial inhomogeneities 
persist even in the optically thin uppermost regions, where the photon bubble instability 
cannot act. Figs. [16] and [T7] show that st rong density inhomoKen e ities a re also present in 
the gas pressure dominated simulation of iHirose. Krolik. fc Stond (120061 ). Photon bubbles 
are almost certainly irrelevant in this case, and therefore cannot be responsible for the 
inhomogeneities that are seen. 



As noted by lHirose. Krolik. fc Stond (120061 ). a more likely explanation for the inhomo- 
geneities is that they are a result of the fact that magnetic forces dominate gas and radiation 
pressure in these layers. Because density only affects gas pressure gradients, and these gen- 
erally play only a small role in the dynamics in these regions, it is not surprising that density 
inhomogeneities develop in response to fluctuating magnetic forces. 



6. Implications for Emergent and Reflected Photons 



Our earlier result that magnetic support extends the upper a tmospheres of accretion 
disks when gas pressure dominates (IHirose. Krolik. fc Stond 120061 ) now appears to apply 
to conditions in which radiation pressure is comparable to gas pressure. This additional 
source of support against gravity creates a lower density photosphere than would otherwise 
be found. As a result, electron scattering is enhanced over absorption, and the matter's 
ionization balance shifts toward more highly- ionized species, weakening absorption edges 
and strengthening emiss ion edges. Both of these effects generically lead to a hardening of 
the emergent spectrum (IBlaes et al.ll2006l ). On the basis of the calculations reported here, 
we expect that magnetically supported, extended atmospheres should continue to exist in 
rings even where radiation and gas pressure are comparable. 

However, we also see strong, irregular density inhomogeneities in the surface layers be- 
neath the scattering photosphere, and these inhomogeneities were also present in the gas 
pressure dominated simulation (Figs. (HHITj). The effect of strong inhomogeneities on the 
emergent spectrum has been investigated in the ordered shock tr ain geometry produced in the 
initia l nonlinear development of the photon bubble instability (IDavis et al.ll2004l ; iBegelman 
20061 ). At least in this geometry, and at least provided thermal absorption/emission dom- 
inates Compton scattering in the energy exchange between radiation and matter, the high 
density regions enhance the thermalization of radiation with matter. This typically produces 
a softening of the spectrum compared to what would emerge from an atmosphere which is a 
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horizontal average of the medium. It remains to be seen which of these two opposing effects 
on the spectrum, softening by inhomogeneities vs. hardening by magnetic support, will be 
dominant, and how much effect Compton scattering can have in those regions where the gas 
temperature rises well above the radiation temperature. 

The very inhomogeneous and irregular nature of the photosphere is also likely to reduce 
the degree of polarization of the emitted photons from a scattering dominated atmosphere, 
as it breaks the plane parallel symmetry. These inhomogenei ties represent a phy s ical r eal- 
ization of the "rough surface" proposed and investigated by IColeman fc Shields! (119901 ) in 
the context of cont i nuum polarization from active galactic nuclei. As first suggested by 
Gnedin &: Silant'evl (119781 ). Faraday depolarization may also be significant. The estimated 
rotation angle at the peak wavelength of a thermal spectrum [i.e., the peak of B\(T)] is 
~ OSttR 1 ^ 2 radians. Here is the Thomson optical depth to the surface where the under- 
lying polarization is imprinted; we expect that because the opacity in these disk atmospheres 
is primarily Thomson scattering, tt — 1. The Faraday rotation is also proportional to R, 
the ratio of magnetic pressure to radiation pressure, if we suppose that the component of the 
field along the line of sight is a fixed and substantial fraction of the total. Using horizontal 
averages of the pressures, the highest total energy epoch in the current simulation (t = 90 
orbits) has R ~ 7 at the lower photosphere and R ~ 3.5 at the upper photosphere. In both, 
the rotation angle is then ~ 2 radians per unit Thomson depth, large enough to give rise to 
significant Faraday depolarization. 

Magnetic support of the surface layers could also have implications for X-ray reflec- 
tion models that incorporate only gas and ra diation pressure in the hydrostatic equilibrium 



[e.g. iNayakshin. Kazanas. Sz Kallmanl 120001 ) . Most notably, the reduced density that fol- 
lows from magnetic support will make photoionization relatively stronger. It is also unclear 
that the discontinuous, thermally stable structures that form in these models would exist if 
magnetic fields dominate the hydrostatic support. The right hand panel of Fig. [TH] shows 
that large density inhomogeneities exist on the scale of a Thomson depth, and these may 
als o affect the reflection spectr u m. T his e ffect has been explored in a prelimin ary fashion 
by iBallantyne. Turner &: Blaed (120041 ) and iBallantyne. Turner. &: Youngl (120051 ). who per- 
formed radiative transfer calculations through one dimensional inhomogeneous structures. 
Density variations of on ly a factor of a few were enough to produce significant changes in 
the reflection spectrum (IBallantyne. Turner fc Blaed 120041 ). and we have more than that at 
the Rosseland mean (~ Thomson) photosphere. Of course, the inhomogeneities themselves 
might be affected by heating by external radiation fields. This topic is clearly one that 
deserves further investigation. 
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7. Summary 

Hydrostatic equilibrium is often thought of in terms of a balance between gravity and 
any of three sorts of pressures: gas, radiation, and magnetic. However, we find that in 
the upper, magnetically-dominated layers of this disk segment, magnetic tension forces are 
almost as strong as magnetic pressure forces, at least in a volume- averaged sense. This is 
clearly because the nonlinear development of the transverse Parker instability is controlling 
the dynamics of the magnetic field in these upper layers. In a mass-averaged sense, however, 
magnetic pressure is still the dominant source of vertical support. 

In contrast to the Parker instability, there is very little evidence for photon bubbles, 
even though their predicted linear growth rates are almost as large as the Parker growth 
rates, and the relevant unstable wavelengths are resolved. 

Strong density inhomogeneities are present in the region between the effective (ther- 
malization) and Rosseland mean (~ scattering) photosph e res. T his was also true in the gas 



pressure dominated simulation of lHirose. Krolik. fc Stond (120061 ). and these inhomogeneities 



are almost certainly the result of the large magnetic forces in the surface layers. Provided 
Compton scattering is less important than true absorption processes in thermally coupling 
photons with matter, these inhomogeneities are likely to enhance thermalization and there- 
fore soften the emergent spectrum. This spectral softening mechanism may counteract to 
some extent the reduction in thermalization, and therefore the spectral hardening, caused by 
the lower average densities that are produced by magnetic support of the atmosphere. Inho- 
mogeneities are likely to reduce the polarization of the emergent photons, and the magnetic 
fields in the atmosphere will reduce it still further through Faraday depolarization. Finally, 
the inhomogeneities we see at the scattering photosphere appear to be strong enough to 
affect X-ray reflection spectra. 

More work needs to be done to elaborate on these conclusions. Given that the size of 
the linearly unstable Parker wavelengths is comparable to the azimuthal size of the box, it 
would be worthwhile in future to run simulations with larger azimuthal domains to see if 
this alters the magnetic equilibrium in the upper layers. Unfortunately, these simulations 
are expensive, so performing a full parameter survey is difficult. 

It is not completely clear why we have failed to find evidence for photon bubbles in the 
present simulation. The accessible growing wavelengths are not too far below the turnover 
wavelength, and increasing the dynamic range of the simulation by decreasing the cell size 
would enable us to access even faster growing wavelengths which might then be better able to 
compete with the Parker instability. However, to fully understand which instability should 
dominate requires, at the very least, a linear instability analysis of an equilibrium with both 
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magnetic pressure support and a diffusive radiation flux. Such an equilibrium would then 
be more representative of the conditions we are finding in the simulation, and would be 
vulnerable to both the Parker and photon bubble instabilities. It could be that these two 
instabilities couple together in nontrivial ways, although their physical driving mechanisms 
are quite different. The transverse Parker instability turns a slow magnetosonic wave into 
an exponentially growing disturbance if magnetic pressure gradients are large enough. The 
fastest growing photon bubble instability is an overstability of, again, a slow magnetosonic 
wave, but is driven by the background radiation pressure. 

Whether the emergent spectrum is harder due to magnetic support or softer due to 
density inhomogeneities when compared to plane parallel atmosphere models that neglect 
magnetic support will require a three dimensional radiative transfer calculation through the 
structures see n in the simulatio n. Monte Carlo calculations may be the way to approach this 



probl em fe.g. iDavis et aJ.ll2Q04f) . and these could also include polarization with Faraday rota- 



tion (lAgol fc Blaeslll996l ). It might also be possible to calculate X-ray reflection in a similar 



fashion, though one would have to include the resulting photoionization self-consistently. 

It remains to be seen what will happen when radiation pressure becomes dominant at 
the midplane. Photon bubbles should have even stronger growth rates in that regime, and 
may be more relevant. However, they will also be more difficult to resolve in a simulation 
because their characteristic wavelengths are only of order the gas pressure scale height. We 
expect to be able to report on simulations with large radiation pressure in the near future. 
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Grants AST-0307657 (OB) and AST-0507455 (JHK). The numerical simulations were carried 
out on the SX8 at the Yukawa Institute for Theoretical Physics of Kyoto University and 
the VPP5000 at the Center for Computational Astrophysics of the National Astronomical 
Observatory of Japan. 



A. Appendix: Simplified Derivation of Parker Instability in a Radiating 

Medium 

Assuming the gas and radiation are locked to the same temperature and that the medium 
is optically thick, the equations of radiation magnetohydro dynamics are 

H + V • (pv) = 0, (Al) 
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P ( % + v ■ Vv ) = " V % + Pg + ■ VB - -L VB 2 + ^F, (A2) 
-|(e + E) + v ■ V(e + £) + f 7 e + ^EJ V • v = -V ■ F, (A3) 

( ' VE, (A4) 



and 



where 



3kfP 

— = V x (v x B), (A5) 

e = . ^ g — r, p s = ^— ^ — , and E = aT A . (A6) 
(7-1) P 

We adopt a simplified, static equilibrium structure in which the gravitational accelera- 
tion g = — gz is constant. Including the linear increase of g with height merely complicates 
the form of the unstable eigenfunctions and is not really justified here given the additional 
simplifications we are making. These include assuming that the background medium is 
isothermal and that the background magnetic energy density is proportional to the density. 
These assumptions imply that the isothermal sound speed in the gas c\ and the Alfven speed 
v a are constants, which greatly simplifies the analysis. A more serious consequence of these 
assumptions is that the background radiation flux F vanishes, and does not help support 
the medium against gravity. Because of this, we immediately lose the driving that causes 
photon bubbles. Nevertheless, we still hope to capture the basic growth rates and wavelength 
scales of the Parker instability in the magnetically dominated upper layers of our simulation. 
Finally, we assume that the background magnetic field is horizontal: B = B(z)y. 

With all these assumptions, hydrostatic equilibrium requires that the background den- 
sity satisfies 

p = p exj) (-Jj:) , (A7) 

where 

H= 2 A±A. (A8 ) 

Linearizing equations (IA3I) and flA4j) about our background, and assuming that our 
disturbances are at short enough wavelengths that the rapid diffusion limit of equation (j3J) 
is satisfied immediately gives us 5T = 0. Assuming the perturbations have a spacetime 
dependence of f(z)exp[i(k x x + k y y — ut)}, and eliminating the magnetic field and pressure 
perturbations using the flux-freezing equation and equation of state, the linearized continuity 
and momentum equations become 

dp r c $ v z 98v z 

- iuj h ik x bv x + ikydVy — + — — = 0, (A9) 

p H oz 



-18- 
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dz 
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(All) 



and 



LO 



l 2 2 

Va + 2#2 



8v, 



3v a cMu* 9 



dz 2 



■ 2d 
2H p HUJCi d~z 



dp 
P 



%k x v\ 
H 



5v x +ik x v\ 



d5v 7 



dz 
(A12) 



We now adopt the ansatz that the z-dependence of the velocity perturbations and rel- 
ative density perturbation is exp[ik z z + z/(2H)], where k z is a real number representing 
the wave vector in the vertical direction. The exp[z/(2H)} factor implies that p|<5w 2 | is con- 
stant with h eight, thereby guaranteein g energy conservation for propagating waves (see e.g. 
section 53 of iMihalas &: Mihalasl Il984 ) . We then obtain the dispersion relation 

kl 



LO 



k 2 + 



AH 2 



(c 2 + v\) + ky A 



+k 2 A 



LO 



AH 2 



+ v\ (2c 2 + vl)(klk 2 




AH 2 



to 



(A13) 



where k 2 = {k 2 + k 2 + k 2 ) 1 / 2 is the total wavenumber of the pertur bation. Equation (1A13|) is 
identical to the original dispersion relation (1.19) of iParkerl (119671 ). restricted to isothermal 
perturbations. 

Instability only exists if k y ^ 0, and when it does, maximum growth occurs for horizontal 
wavenumbers ik z = 0). In this case, instability occurs only for wavelengths X y = 2it/k y 
satisfying 

2irci(2 C ?+vp . f . _n 

9 ( C 2+l,2)l/2 11 I^X - U, 

— 1 — it k r — > OO. 



V ^* A Parker 



(A14) 



Note that if the magnetic energy density is dominant over gas pressure, Ap ar k cr — ^c{v^jg. 
The wavelength X y for maximum growth occurs at 



2| 1/2 



A, 



2^(2cf+v 2 A )\v 2 A -cf\c[ 

9{(c?+^i)[2 1 / 2 (c 2 +i-i) 3 /2- Ci ( C 2+3t;2)]}i/2 
2tt»;2 ( 2c 2 + „2)i/2 e i/2 



if kr = 0, 



if k n 



(A15) 



00. 



9[2 1 /2(2 C 2 + ^ ) 3/2_(4 c 2 + 3^ )c .]l/2 

In the limit of strong magnetic fields, A max ~ 2nv ^v^Ci) 1 ! 2 / (2 l l 4 g). The corresponding 
maximum growth rate is 
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For strong magnetic fields, this becomes simply 7 max ~ g/v\. 
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Fig. 1. — Wavelength Ar below which rapid photon diffusion applies (lower of each pair of 
curves), and wavelength \ s above which photons are trapped (upper of each pair of curves) 
for acoustic disturbances in the simulation. The solid and dashed pair of curves correspond 
to horizontally averaged conditions at the highest (t = 90 orbits) and lowest (t = 150 
orbits) total energy epochs, respectively. The curves extend only out to the Rosseland mean 
photosphere of the horizontally averaged structure in each case, outside of which photon 
diffusion does not apply. For comparison, the lower horizontal dot-dashed line shows the 
length 8r A0 of eight grid zones in the azimuthal direction, close to the minimum wavelength 
that the code can resolve. The upper horizontal dot-dashed line shows the length 3H of the 
box in the azimuthal direction. 
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Fig. 2. — Scaled photon/gas thermal coupling frequency calculated from the horizontally 
averaged structure at t = 90 orbits (solid) and t = 150 orbits (dashed). Tight thermal 
coupling is a good approximation for photon bubbles if this scaled frequency is much greater 
than unity. 
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Fig. 3. — Range of wavelengths over which the asymptotic photon bubble growth rates of 
Fig. H] should be valid (provided these wavelengths are also in the rapid photon diffusion 
regime), as a function of height in the box at t = 90 orbits (solid) and t = 150 orbits 
(dashed). The solid and upper dashed curves in each epoch correspond to the turnover 
wavelength, while the lower dashed curve on the right corresponds to the thermal cutoff 
wavelength at t = 150 orbits. (The thermal cutoff wavelength for t = 90 orbits is less than 
0.1H at all heights, and is also less than 0.1H for z < 2.8H at t = 150 orbits.) Rapid 
growth rates should exist between the thermal cutoff and turnover wavelengths. The lower 
and upper horizontal dashed lines again indicate eight cell sizes 8r and the length 3H of 
the box in the azimuthal direction, respectively. 
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Fig. 4. — Short wavelength photon bubble growth rate 7^ scaled with the orbital period 
P or b = 2-k/VL, as a function of height at t = 90 orbits (solid) and t = 150 orbits (dashed). 
Growth rates are shown only at heights deeper than the photosphere of the horizontally 
averaged structure at that epoch. 
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Fig. 5. — Maximum Parker instability growth rate 7 max from equation flA16j) . scaled with 
the orbital period P or b = 27r/Q, as a function of height in the box at t = 90 orbits (solid) 
and t = 150 orbits (dashed). The two nearly identical curves for each epoch correspond to 
the k x = and oo limits. 
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Fig. 6. — Minimum Parker instability wavelength from equation (IA14I) (lower of each pair 
of curves) and wavelength for maximal Parker growth rate from equation (IA15I) (upper of 
each pair of curves) as a function of height in the box at t = 90 orbits (solid curve pair) and 
t = 150 orbits (dashed curve pair). Eight cell sizes and the length of the box in the azimuthal 
direction are indicated again by the lower and upper horizontal dashed lines, respectively. 
The left hand panel corresponds to k x = 0, while the right hand panel corresponds to 
k x — > oo. The dependence on k x is very weak. 



-27- 




Fig. 7. — Density in regions under the horizontally-averaged upper photosphere (dashed 
line) in a fixed radial slice near the middle of the box at t = 90 orbits. The arrows show the 
projections of magnetic field vectors (black) and velocity vectors relative to the background 
shear flow (red) in this radial slice, computed at the position of the tail of each arrow. 
Horizontal dashed lines indicate the positions of the Rosseland mean photosphere (upper) 
and effective photosphere (lower), computed from the horizontally averaged structure. 
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Fig. 8. — Distribution of cosines of the angle 9 between the perturbed radiation flux vector 
and the fluid velocity relative to the background shear flow, over horizontal slices in the 
simulation at t = 90 orbits. Different histograms correspond to different heights: +2AH 
(black), +2.9H (red), +3.3H (green), and +3.8H (blue). 
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Fig. 9. — Spatially horizontally averaged (left) and mass-weighted horizontally averaged 
(right) contributions to the outward vertical accelerations from magnetic pressure gradients 
(dashed) and magnetic tension (dotted), at t = 90 orbits. The solid curves show the sum of 
these magnetic accelerations. 




Fig. 10. — Vertical magnetic tension force per unit volume vs. density at each grid zone in 
a, z = 2.75H horizontal slice through the simulation at t = 90 orbits (left) and in a z = 2H 
horizontal slice through the simulation at t — 150 orbits (right). The tension has been scaled 
by the average magnitude of the tension force in this horizontal slice, and the density has 
been scaled by the average density. 
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Fig. 11. — Density in the upper layers (2H < z < 6H) in a fixed radial slice near the 
middle of the box at time 150 orbits. The arrows show the projections of magnetic field 
vectors (black) and velocity vectors relative to the background shear flow (red) into this 
radial slice, computed at the position of the tail of each arrow. The Rosseland mean and 
effective photospheres of the horizontally averaged structure at this time are indicated by 
the upper and lower dashed lines, respectively. 
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Fig. 12. — Magnetic field lines in the simulation domain at 90 orbits (left) and 150 orbits 
(right). The lines are color coded with the local fluid density scaled by the horizontally 
averaged density at the same height. 
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Fig. 13. — The fractional horizontal area over which higher (lower) than average densities 
are associated with vertically inward (outward) magnetic tension forces, as a function of 
height in the simulation. The left panel is for t = 90 orbits and the right panel is for t = 150 
orbits. 
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Fig. 14. — Density in the uppermost surface layers (3H < z < 8H) in a fixed radial slice 
near the middle of the box at time 60 orbits from the gas pressure dominated simulation of 



Hirose. Krolik. &: Stone! (120061 ) . The arrows show the projections of magnetic field vectors 



(black) and velocity vectors relative to the background shear flow (red) into this radial 
slice, computed at the position of the tail of each arrow. The Rosseland mean and effective 
photospheres of the horizontally averaged structure at this time are indicated by the upper 
and lower dashed lines, respectively. 
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Fig. 16. — Distribution of densities, scaled with the horizontally averaged density, at the 
heights of the upper effective (left) and Rosseland mean (right) photospheres at 90 orbits 
(solid histograms) and 150 orbits (dashed histograms). For comparison, the dotted his- 
tograms show the corresponding densit y distributions in the gas pressure dominated simu- 
lation of iHirose. Krolik. fc Stone J2006h at 60 orbits. 
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Fig. 17. — Root mean square fractional fluctuation of density as a function of height at 90 
orbits (solid) and 150 orbits (dashed). For co mparison, the dotted curve sho ws the same 
thing in the gas pressure dominated simulation of lHirose. Krolik. & Stone! (120061 ) at 60 orbits. 
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Fig. 18. — Ratio of gas to radiation temperature near the upper photosphere in a fixed 
radial slice near the middle of the box at t=90 orbits. (The radial slice is identical to 
that chosen in Fig. [7J) Horizontal dashed lines indicate the positions of the Rosseland 
mean photosphere (upper) and effective photosphere (lower), computed from the horizontally 
averaged structure. 
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Fig. 19. — Same as Fig. [18] except at t=150 orbits. (The corresponding density structure is 
shown in Fig. [TT1 ) 



